library(AlgDesign)
library(plyr)
library(Matrix)
library(mvtnorm)
library(mlogit)
#################################
# Step 1: Simulate data
#################################

ni <- 300 # num of individuals
np <- 5   # num of purchases for each individual
n <- ni * np
K <- 30  # num of brands
p <- 2   # number of individual specific vars
nG <- 7
nL <- 3
rho <- 0.3

sigma <- matrix(0, nG, nG)
sigma <- rho^abs(col(sigma) - row(sigma))

set.seed(11)
# data structure
sid <- rep(1:ni, rep(np, ni))
id <- rep(1:n, rep(K, n))
da <- data.frame(id = id, 
                 sid =  rep(sid, rep(K, n)),
                 choice = rep(c(1, rep(0, K - 1)), n),
                 brand = rep(1:K, n))   

#brand specific vars
ss <- matrix(0, 5, 5)
ss <- rho^abs(col(ss) - row(ss))
x <- rmvnorm(n * K, sigma = ss)
x <- as.data.frame(x)
names(x) <- paste("x", 1:5, sep = "")
da <- data.frame(da, x)

# individual specific variables
da2 <- matrix(rnorm(ni * p), ni, p)
da2 <- data.frame(sid = 1:ni, da2)
names(da2)[-1] <- paste0("x", 6:(5 + p))

#add in individual specific vars
da <- join(da, da2) 

# generate attributes
X <- rmvnorm(K, sigma = sigma)
X2 <- apply(X, 2, function(x){
  xx <- cut(x, nL)
  levels(xx) <- 1:nL
  return(xx)
})
X2 <- as.data.frame(X2)
# check if colinear
tmp <- model.matrix(~., as.data.frame(X2))
if (qr(tmp)$rank != ncol(tmp))
  stop("X does not have column rank")

Z <- do.call(rbind, rlply(ni * np, X2))
names(Z) <- paste("z", 1:nG, sep = "")
da <- cbind(da, data.frame(Z))

da <- transform(da, id = factor(id), 
                brand = factor(brand),
                sid = factor(sid))


X <- model.matrix(~ z1 + z2 + z3 + z4 + z5 + x1 + x2, da)[, -1]
if (qr(X)$rank != ncol(X))
  stop("X does not have column rank")

# convert for object "mlogit.data"
da <- mlogit.data(da, choice = "choice", shape = "long", alt.var = "brand", id.var = "sid")